library(mixOmics)
library(ggplot2)

##读入文件
#门水平丰度表
phylum <- read.delim('phylum_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
phylum <- data.frame(t(phylum))

#样本分组文件
group <- read.delim('group.txt', sep = '\t', stringsAsFactors = FALSE)

##PLS-DA 分析
#基于门水平丰度表，只展示前 3 个排序轴
phylum <- phylum[group$names, ]

head(phylum)
head(group$group)
plsdaRs <-plsda(phylum, group$group, ncomp = 3)

#简要查看结果
plsda_result
#或
names(plsda_result)

#查看排序轴解释量
plsda_result$explained_variance$X
#查看样本排序坐标
plsda_result$variates$X
#查看细菌门类群排序坐标
plsda_result$loadings$X

#使用 plotIndiv() 绘制 PLS-DA 分析结果
plotIndiv(plsda_result, ind.names = TRUE, style = 'ggplot2')
plotIndiv(plsda_result, ind.names = TRUE, style = '3d')

#提取坐标轴解释量（前两轴）
plsda_result_eig <- {plsda_result$explained_variance$X}[1:2]

#提取样本点坐标（前两轴）
sample_site <- data.frame(plsda_result$variates)[1:2]

#为样本点坐标添加分组信息
sample_site$names <- rownames(sample_site)
names(sample_site)[1:2] <- c('plsda1', 'plsda2')
sample_site <- merge(sample_site, group, by = 'names', all.x = TRUE)

#可选输出各样本的 PLS-DA 分析结果
write.table(sample_site, 'plsda_sample.txt', row.names = FALSE, sep = '\t', quote = FALSE)

#使用 ggplot2 简单绘制 PLS-DA 结果图
plsda_plot <- ggplot(sample_site, aes(plsda1, plsda2, color = group, label = names)) +
    geom_point(size = 1.5, alpha = 0.6) +
    stat_ellipse(show.legend = FALSE) +	#添加 95% 置信椭圆
    scale_color_manual(values = c('#1D7ACC', '#F67433', '#00815F')) +
    theme(panel.grid = element_line(color = 'grey50'), panel.background = element_rect(color = 'black', fill = 'transparent')) +
    theme(legend.title = element_blank(), legend.key = element_rect(fill = 'transparent')) +
    labs(x = paste('PLS-DA axis1 ( explained variance ', round(100 * plsda_result_eig[1], 2), '% )', sep = ''), y = paste('PLS-DA axis2 ( explained variance ', round(100 * plsda_result_eig[2], 2), '% )', sep = ''))

#ggsave(limitsize = FALSE,'plsda_plot.pdf', plsda_plot, width = 6, height = 5)
ggsave(limitsize = FALSE,'plsda_plot.png', plsda_plot, width = 6, height = 5)
